MODULE GSBCOST_GRAD
!*************************************************
! CALCULATE COST FUNCTION AND THE RELEVENT GRADIENTS OF GS-BALANCE PENALTY TERM, AND THE PENALTY TERM OF CONTINUAL CONDITION.
! HISTORY: MODIFIED AND DEPARTED FROM STMAS4D_CORE MODULE BY ZHONGJIE HE, JANUARY 2008.
!*************************************************

  USE PRMTRS_STMAS
  USE GENERALTOOLS, ONLY : GDERIVELB, GDERIVEIT, GDERIVERB

  PUBLIC  GSBLNCOST_NON_UNIFORM_Z, GSBLNGRAD_NON_UNIFORM_Z, GSBLNCOST_P_HS, GSBLNGRAD_P_HS, CNTNSCOST, CNTNSGRAD, GSBLNCOST_P, GSBLNGRAD_P

!***************************************************
!!COMMENT:
!   THIS MODULE IS USED BY THE MODULE OF STMAS4D_CORE TO CALCULATE THE COST FUNCTION AND THE RELEVENT GRADIENTS OF PENALTY TERMS.
!   SUBROUTINES:
!      GSBLNCOST_NON_UNIFORM_Z : CALCULATE THE COSTFUNCTION OF GS-BALANCE FOR THE SIGMA AND HEIGHT COORDINATE SYSTEM.
!      GSBLNGRAD_NON_UNIFORM_Z : CALCULATE GRADIENTS OF COSTFUNCTIONS CORRESPOND TO GSBLNCOST_NON_UNIFORM_Z.
!      GSBLNCOST_P_HS          : CALCULATE THE COSTFUNCTION FOR THE CASE OF PRESSURE COORDINATE, THE HYDROSTATIC APPROXIMATION IS USED.
!      GSBLNGRAD_P_HS          : CALCULATE GRADIENTS OF COSTFUNCTIONS CORRESPOND TO GSBLNCOST_P_HS.
!      CNTNSCOST               : CALCULATE THE PENALTY TERM OF CONTINUAL CONDITIONS.
!      CNTNSGRAD               : CALCULATE THE GRADIENT OF COSTFUNCTIONS CORRESPOND TO CNTNSCOST.
!      GSBLNCOST_P             : JUST LIKE GSBLNCOST_P_HS, NOT USE THE HYDROSTATIC APPROXIMATION.
!      GSBLNGRAD_P             : CALCULATE GRADIENTS OF COSTFUNCTIONS CORRESPOND TO GSBLNCOST_P.
CONTAINS

SUBROUTINE GSBLNCOST_NON_UNIFORM_Z
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF COST FUNCTION
! HISTORY: AUGUST 2007, CODED by WEI LI.
!          FEBRUARY 2008, MODIFIED BY ZHONGJIE HE
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,PP
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN

! CALCULATE THE TERM OF DP/DX=DEN*COR*V
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=ZZZ(I,J,K1,T)
      Z2=ZZZ(I,J,K2,T)
      Z3=ZZZ(I,J,K3,T)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)

      GS=((GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
         /(XXX(I2,J)-XXX(I1,J)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J)))*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
        -((GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
         /(YYY(I,J2)-YYY(I,J1)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1)))*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
         -COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))

      COSTFUN=COSTFUN+PNLT_PV*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF

! CALCULATE THE TERM OF DP/DY=-DEN*COR*U
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=ZZZ(I,J,K1,T)
      Z2=ZZZ(I,J,K2,T)
      Z3=ZZZ(I,J,K3,T)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)

      GS=((GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
         /(YYY(I,J2)-YYY(I,J1)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1)))*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
        +((GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
         /(XXX(I2,J)-XXX(I1,J)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J)))*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
         +COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))

      COSTFUN=COSTFUN+PNLT_PU*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNCOST_NON_UNIFORM_Z

SUBROUTINE GSBLNGRAD_NON_UNIFORM_Z
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF GRADIENT
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,PP
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=ZZZ(I,J,K1,T)
      Z2=ZZZ(I,J,K2,T)
      Z3=ZZZ(I,J,K3,T)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)

      GS=((GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
         /(XXX(I2,J)-XXX(I1,J)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J)))*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
        -((GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
         /(YYY(I,J2)-YYY(I,J1)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1)))*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
         -COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))
!
      GRADINT(I,J,K,T,VV) =GRADINT(I,J,K,T,VV) +2.0*PNLT_PV*GS*COR(I,J)*DEN(I,J,K,T)*(-1.0)
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         /(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         /(XXX(I2,J)-XXX(I1,J))*(-1.0)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *AZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *BZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *CZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)-2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)-2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1))*(-1.0)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)-2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *AZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)-2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *BZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)-2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *CZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))

    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=ZZZ(I,J,K1,T)
      Z2=ZZZ(I,J,K2,T)
      Z3=ZZZ(I,J,K3,T)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      GS=((GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
         /(YYY(I,J2)-YYY(I,J1)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1)))*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
        +((GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
         /(XXX(I2,J)-XXX(I1,J)) &
         -(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
         *(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J)))*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
         +COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))
!
      GRADINT(I,J,K,T,UU) =GRADINT(I,J,K,T,UU) +2.0*PNLT_PU*GS*COR(I,J)*DEN(I,J,K,T)
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)+2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)+2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1))*(-1.0)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *AZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *BZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PU*GS*SCL(PP)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *CZ*(-1.0)*(ZZZ(I,J2,K,T)-ZZZ(I,J1,K,T))/(YYY(I,J2)-YYY(I,J1))
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         /(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         /(XXX(I2,J)-XXX(I1,J))*(-1.0)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *AZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *BZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PV*GS*SCL(PP)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *CZ*(-1.0)*(ZZZ(I2,J,K,T)-ZZZ(I1,J,K,T))/(XXX(I2,J)-XXX(I1,J))
!
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNGRAD_NON_UNIFORM_Z


SUBROUTINE GSBLNCOST_P_HS
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF COST FUNCTION
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  REAL ,PARAMETER :: GA=9.8D0
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,UU,VV,PP
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN

! CALCULATE THE TERM OF DP/DX=DEN*COR*V
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      GS=(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(XXX(I2,J)-XXX(I1,J))*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(VV)*GA*1.0 &
        -(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(YYY(I,J2)-YYY(I,J1))*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(VV)*GA*0.0 &
        -COR(I,J)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))
      COSTFUN=COSTFUN+PNLT_PV*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDIF

! CALCULATE THE TERM OF DP/DY=-DEN*COR*U
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      GS=(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(YYY(I,J2)-YYY(I,J1))*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(UU)*GA*1.0 &
        +(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(XXX(I2,J)-XXX(I1,J))*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(UU)*GA*0.0 &
        +COR(I,J)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))
      COSTFUN=COSTFUN+PNLT_PU*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNCOST_P_HS

SUBROUTINE GSBLNGRAD_P_HS
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF GRADIENT
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  REAL ,PARAMETER :: GA=9.8D0
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,UU,VV,PP
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      GS=(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(XXX(I2,J)-XXX(I1,J))*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(VV)*GA*1.0 &
        -(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(YYY(I,J2)-YYY(I,J1))*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(VV)*GA*0.0 &
        -COR(I,J)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))
!
      GRADINT(I,J,K,T,VV) =GRADINT(I,J,K,T,VV) +2.0*PNLT_PV*GS*COR(I,J)*(-1.0)
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(VV)*1.0 &
	                     /(XXX(I2,J)-XXX(I1,J))*GA
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(VV)*(-1.0)*1.0 &
                         /(XXX(I2,J)-XXX(I1,J))*GA
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)-2.0*PNLT_PU*GS*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(VV)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1))*GA
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)-2.0*PNLT_PU*GS*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(VV)*(-1.0)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1))*GA
!
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      GS=(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(YYY(I,J2)-YYY(I,J1))*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(UU)*GA*1.0 &
        +(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(XXX(I2,J)-XXX(I1,J))*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(UU)*GA*0.0 &
        +COR(I,J)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))
!
      GRADINT(I,J,K,T,UU) =GRADINT(I,J,K,T,UU) +2.0*PNLT_PU*GS*COR(I,J)
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)+2.0*PNLT_PU*GS*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(UU)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1))*GA
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)+2.0*PNLT_PU*GS*(SCL(PP)*Z_TRANS)/SCP(YSL)/SCP(CSL)/SCL(UU)*(-1.0)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1))*GA
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(UU)*0.0 &
	                     /(XXX(I2,J)-XXX(I1,J))*GA
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*(SCL(PP)*Z_TRANS)/SCP(XSL)/SCP(CSL)/SCL(UU)*(-1.0)*0.0 &
                         /(XXX(I2,J)-XXX(I1,J))*GA
!
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNGRAD_P_HS

SUBROUTINE GSBLNCOST_P
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF COST FUNCTION
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,PP
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      GS=(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(XXX(I2,J)-XXX(I1,J))*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
        -(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(YYY(I,J2)-YYY(I,J1))*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
        +COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))
      COSTFUN=COSTFUN+PNLT_PV*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDIF
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      GS=(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(YYY(I,J2)-YYY(I,J1))*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
        +(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(XXX(I2,J)-XXX(I1,J))*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
        -COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))
      COSTFUN=COSTFUN+PNLT_PU*GS**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNCOST_P

SUBROUTINE GSBLNGRAD_P
!*************************************************
! GEOSTROPHIC BALANCE PENALTY TERM OF GRADIENT
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,PP
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: GS
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  PP=PRESSURE
  IF(PNLT0PU.LT.1.0E-10.AND.PNLT0PV.LT.1.0E-10)RETURN
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      GS=(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(XXX(I2,J)-XXX(I1,J))*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
        -(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(YYY(I,J2)-YYY(I,J1))*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
        +COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,VV)+GRDBKGND(I,J,K,T,VV))
!
      GRADINT(I,J,K,T,VV) =GRADINT(I,J,K,T,VV) +2.0*PNLT_PV*GS*COR(I,J)*DEN(I,J,K,T)
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                             /(XXX(I2,J)-XXX(I1,J)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*(-1.0)*1.0 &
                         /(XXX(I2,J)-XXX(I1,J)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*AZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*BZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*1.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*CZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)-2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)-2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*(-1.0)*0.0 &
                         /(YYY(I,J2)-YYY(I,J1)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)-2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*AZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)-2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*BZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)-2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(VV)/SCP(DSL)*0.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*CZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      GS=(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP)+GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(YYY(I,J2)-YYY(I,J1))*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
        +(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP)+GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP)) &
        /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
         +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ) &
        /(XXX(I2,J)-XXX(I1,J))*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
        -COR(I,J)*DEN(I,J,K,T)*(GRDANALS(I,J,K,T,UU)+GRDBKGND(I,J,K,T,UU))
!
      GRADINT(I,J,K,T,UU) =GRADINT(I,J,K,T,UU) +2.0*PNLT_PU*GS*COR(I,J)*DEN(I,J,K,T)*(-1.0)
!
      GRADINT(I,J2,K,T,PP)=GRADINT(I,J2,K,T,PP)+2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J1,K,T,PP)=GRADINT(I,J1,K,T,PP)+2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*(-1.0)*1.0 &
                         /(YYY(I,J2)-YYY(I,J1)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*AZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*BZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PU*GS*SCP(PSL)/SCP(YSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*1.0 &
                         *(GRDANALS(I,J2,K,T,PP)-GRDANALS(I,J1,K,T,PP) &
                          +GRDBKGND(I,J2,K,T,PP)-GRDBKGND(I,J1,K,T,PP))/(YYY(I,J2)-YYY(I,J1))*CZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I2,J,K,T,PP)=GRADINT(I2,J,K,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                             /(XXX(I2,J)-XXX(I1,J)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I1,J,K,T,PP)=GRADINT(I1,J,K,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*(-1.0)*0.0 &
                         /(XXX(I2,J)-XXX(I1,J)) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)
!
      GRADINT(I,J,K1,T,PP)=GRADINT(I,J,K1,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*AZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K2,T,PP)=GRADINT(I,J,K2,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*BZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
      GRADINT(I,J,K3,T,PP)=GRADINT(I,J,K3,T,PP)+2.0*PNLT_PV*GS*SCP(PSL)/SCP(XSL)/SCP(CSL)/SCL(UU)/SCP(DSL)*0.0 &
                         *(GRDANALS(I2,J,K,T,PP)-GRDANALS(I1,J,K,T,PP) &
                          +GRDBKGND(I2,J,K,T,PP)-GRDBKGND(I1,J,K,T,PP))/(XXX(I2,J)-XXX(I1,J))*CZ*(-1.0) &
                         /(GRDANALS(I,J,K1,T,PP)*AZ+GRDANALS(I,J,K2,T,PP)*BZ+GRDANALS(I,J,K3,T,PP)*CZ &
                          +GRDBKGND(I,J,K1,T,PP)*AZ+GRDBKGND(I,J,K2,T,PP)*BZ+GRDBKGND(I,J,K3,T,PP)*CZ)**2
!
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE GSBLNGRAD_P


SUBROUTINE CNTNSCOST
!*************************************************
! CONTINOUS EQUATION
! HISTORY: OCTOBER 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,WW
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: CE
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  WW=W_CMPNNT
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      CE=(GRDANALS(I2,J,K,T,UU)-GRDANALS(I1,J,K,T,UU)+GRDBKGND(I2,J,K,T,UU)-GRDBKGND(I1,J,K,T,UU)) &
        /(XXX(I2,J)-XXX(I1,J))*SCL(UU)/SCP(XSL)*SCP(PSL)/SCL(WW)*1.0 &
        -(GRDANALS(I,J2,K,T,UU)-GRDANALS(I,J1,K,T,UU)+GRDBKGND(I,J2,K,T,UU)-GRDBKGND(I,J1,K,T,UU)) &
        /(YYY(I,J2)-YYY(I,J1))*SCL(UU)/SCP(YSL)*SCP(PSL)/SCL(WW)*0.0 &
        +(GRDANALS(I,J2,K,T,VV)-GRDANALS(I,J1,K,T,VV)+GRDBKGND(I,J2,K,T,VV)-GRDBKGND(I,J1,K,T,VV)) &
        /(YYY(I,J2)-YYY(I,J1))*SCL(VV)/SCP(YSL)*SCP(PSL)/SCL(WW)*1.0 &
        +(GRDANALS(I2,J,K,T,VV)-GRDANALS(I1,J,K,T,VV)+GRDBKGND(I2,J,K,T,VV)-GRDBKGND(I1,J,K,T,VV)) &
        /(XXX(I2,J)-XXX(I1,J))*SCL(VV)/SCP(XSL)*SCP(PSL)/SCL(WW)*0.0 &
        +(GRDANALS(I,J,K1,T,WW)*AZ+GRDANALS(I,J,K2,T,WW)*BZ+GRDANALS(I,J,K3,T,WW)*CZ &
         +GRDBKGND(I,J,K1,T,WW)*AZ+GRDBKGND(I,J,K2,T,WW)*BZ+GRDBKGND(I,J,K3,T,WW)*CZ)
      COSTFUN=COSTFUN+PNLT0PU*CE**2
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE CNTNSCOST

SUBROUTINE CNTNSGRAD
!*************************************************
! CONTINOUS EQUATION
! HISTORY: OCTOBER 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,I1,I2,J1,J2,K1,K2,K3,UU,VV,WW
  REAL     :: Z1,Z2,Z3,AZ,BZ,CZ
  REAL     :: CE
! --------------------

  UU=U_CMPNNT
  VV=V_CMPNNT
  WW=W_CMPNNT
  IF(NUMGRID(1).GE.2.AND.NUMGRID(2).GE.2.AND.NUMGRID(3).GE.3.AND.NUMSTAT.GE.3)THEN
    DO T=1,NUMGRID(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      I1=MAX0(I-1,1)
      I2=MIN0(I+1,NUMGRID(1))
      J1=MAX0(J-1,1)
      J2=MIN0(J+1,NUMGRID(2))
      K1=MAX0(K-1 ,1)
      K1=MIN0(K1  ,NUMGRID(3)-2)
      K2=MIN0(K1+1,NUMGRID(3)-1)
      K3=MIN0(K2+1,NUMGRID(3))
      Z1=PPP(K1)
      Z2=PPP(K2)
      Z3=PPP(K3)
      IF(K.EQ.1)CALL GDERIVELB(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.GE.2.AND.K.LE.NUMGRID(3)-1)CALL GDERIVEIT(Z1,Z2,Z3,AZ,BZ,CZ)
      IF(K.EQ.NUMGRID(3))CALL GDERIVERB(Z1,Z2,Z3,AZ,BZ,CZ)
      CE=(GRDANALS(I2,J,K,T,UU)-GRDANALS(I1,J,K,T,UU)+GRDBKGND(I2,J,K,T,UU)-GRDBKGND(I1,J,K,T,UU)) &
        /(XXX(I2,J)-XXX(I1,J))*SCL(UU)/SCP(XSL)*SCP(PSL)/SCL(WW)*1.0 &
        -(GRDANALS(I,J2,K,T,UU)-GRDANALS(I,J1,K,T,UU)+GRDBKGND(I,J2,K,T,UU)-GRDBKGND(I,J1,K,T,UU)) &
        /(YYY(I,J2)-YYY(I,J1))*SCL(UU)/SCP(YSL)*SCP(PSL)/SCL(WW)*0.0 &
        +(GRDANALS(I,J2,K,T,VV)-GRDANALS(I,J1,K,T,VV)+GRDBKGND(I,J2,K,T,VV)-GRDBKGND(I,J1,K,T,VV)) &
        /(YYY(I,J2)-YYY(I,J1))*SCL(VV)/SCP(YSL)*SCP(PSL)/SCL(WW)*1.0 &
        +(GRDANALS(I2,J,K,T,VV)-GRDANALS(I1,J,K,T,VV)+GRDBKGND(I2,J,K,T,VV)-GRDBKGND(I1,J,K,T,VV)) &
        /(XXX(I2,J)-XXX(I1,J))*SCL(VV)/SCP(XSL)*SCP(PSL)/SCL(WW)*0.0 &
        +(GRDANALS(I,J,K1,T,WW)*AZ+GRDANALS(I,J,K2,T,WW)*BZ+GRDANALS(I,J,K3,T,WW)*CZ &
         +GRDBKGND(I,J,K1,T,WW)*AZ+GRDBKGND(I,J,K2,T,WW)*BZ+GRDBKGND(I,J,K3,T,WW)*CZ)
      GRADINT(I2,J,K,T,UU)=GRADINT(I2,J,K,T,UU)+2.0*PNLT0PU*CE &
                         /(XXX(I2,J)-XXX(I1,J))*SCL(UU)/SCP(XSL)*SCP(PSL)/SCL(WW)*1.0
      GRADINT(I1,J,K,T,UU)=GRADINT(I1,J,K,T,UU)+2.0*PNLT0PU*CE*(-1.0) &
                         /(XXX(I2,J)-XXX(I1,J))*SCL(UU)/SCP(XSL)*SCP(PSL)/SCL(WW)*1.0
      GRADINT(I,J2,K,T,UU)=GRADINT(I,J2,K,T,UU)+2.0*PNLT0PU*CE*(-1.0) &
                         /(YYY(I,J2)-YYY(I,J1))*SCL(UU)/SCP(YSL)*SCP(PSL)/SCL(WW)*0.0
      GRADINT(I,J1,K,T,UU)=GRADINT(I,J1,K,T,UU)+2.0*PNLT0PU*CE &
                         /(YYY(I,J2)-YYY(I,J1))*SCL(UU)/SCP(YSL)*SCP(PSL)/SCL(WW)*0.0
      GRADINT(I,J2,K,T,VV)=GRADINT(I,J2,K,T,VV)+2.0*PNLT0PU*CE &
                         /(YYY(I,J2)-YYY(I,J1))*SCL(VV)/SCP(YSL)*SCP(PSL)/SCL(WW)*1.0
      GRADINT(I,J1,K,T,VV)=GRADINT(I,J1,K,T,VV)+2.0*PNLT0PU*CE*(-1.0) &
                         /(YYY(I,J2)-YYY(I,J1))*SCL(VV)/SCP(YSL)*SCP(PSL)/SCL(WW)*1.0
      GRADINT(I2,J,K,T,VV)=GRADINT(I2,J,K,T,VV)+2.0*PNLT0PU*CE &
                         /(XXX(I2,J)-XXX(I1,J))*SCL(VV)/SCP(XSL)*SCP(PSL)/SCL(WW)*0.0
      GRADINT(I1,J,K,T,VV)=GRADINT(I1,J,K,T,VV)+2.0*PNLT0PU*CE*(-1.0) &
                         /(XXX(I2,J)-XXX(I1,J))*SCL(VV)/SCP(XSL)*SCP(PSL)/SCL(WW)*0.0
      GRADINT(I,J,K1,T,WW)=GRADINT(I,J,K1,T,WW)+2.0*PNLT0PU*CE*AZ
      GRADINT(I,J,K2,T,WW)=GRADINT(I,J,K2,T,WW)+2.0*PNLT0PU*CE*BZ
      GRADINT(I,J,K3,T,WW)=GRADINT(I,J,K3,T,WW)+2.0*PNLT0PU*CE*CZ
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  RETURN
END SUBROUTINE CNTNSGRAD

END MODULE GSBCOST_GRAD
